function J = jacobian(func, x, y, a, b, c, d)
% jacobian 计算给定模型函数 func 在状态 (x, y) 处的雅可比矩阵
%
% 输入：
%   func - 模型函数句柄，形如 [x_next, y_next] = func(x, y, a, b, c, d)
%   x, y - 当前状态值
%   a, b, c, d - 模型参数
%
% 输出：
%   J - 2x2 雅可比矩阵，其中
%         J(1,1) = ∂(x_next)/∂x,  J(1,2) = ∂(x_next)/∂y,
%         J(2,1) = ∂(y_next)/∂x,  J(2,2) = ∂(y_next)/∂y

    % 设置微小扰动步长
    h = 1e-6;
    
    % 计算原始函数值
    [f0x, f0y] = func(x, y, a, b, c, d);
    
    % 计算在 x 方向扰动后的函数值
    [fx1, fy1] = func(x + h, y, a, b, c, d);
    
    % 计算在 y 方向扰动后的函数值
    [fx2, fy2] = func(x, y + h, a, b, c, d);
    
    % 有限差分计算偏导数
    dfdx = ([fx1; fy1] - [f0x; f0y]) / h;
    dfdy = ([fx2; fy2] - [f0x; f0y]) / h;
    
    % 组装雅可比矩阵
    J = [dfdx, dfdy];
end
